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Abstract 


Radiation symmetry evaluation is critical to the laser driven Inertial Confinement Fusion (ICF), which is usually done by 
solving a view-factor equation model. The model is nonlinear, and the number of equations can be very large when the size of 
discrete mesh element is very small to achieve a prescribed accuracy, which may lead to an intensive equation solving process. In 
this paper, an efficient radiation symmetry analysis approach based on sparse representation is presented, in which, 1) the Spherical 
harmonics, annular Zernike polynomials and Legendre-Fourier polynomials are employed to sparsely represent the radiation flux 
on the capsule and cylindrical cavity, and the nonlinear energy equilibrium equations are transformed into the equations with 
sparse coefficients, which means there are many redundant equations, 2) only a few equations are selected to recover such sparse 
coefficients with Latin hypercube sampling, 3) a Conjugate Gradient Subspace Thresholding Pursuit (CGSTP) algorithm is then 
given to rapidly obtain such sparse coefficients equation with as few iterations as possible. Finally, the proposed method is validated 
with two experiment targets for Shenguang II and Shenguang III laser facility in China. The results show that only one tenth 
of computation time is required to solve one tenth of equations to achieve the radiation flux with comparable accuracy. Further 
more, the solution is much more efficient as the size of discrete mesh element decreases, in which, only 1.2% computation time 
is required to obtain the accurate result. 
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I. INTRODUCTION 


ONTROLLABLE nuclear fusion is potential to solve the energy crisis in the future, and laser driven Inertial Confinement 

Fusion (ICF) is supposed to be one of the most promising ways [2]. To achieve this goal, the driven asymmetry on 
the centrally located capsule needs to be evaluated and remains no more than 1% during the fusion process. Due to limited 
laser beams for NIF laser facility, i.e., 2 laser entrance holes and a few discrete diagnose holes, the radiation asymmetry 
may be much larger than such prescribed threshold [3]. Therefore, we need to efficiently compute the driven flux reached 
the capsule to evaluate its driven asymmetry. In the ICF, the radiation flux on the fuel capsule is related to the laser-plasma 
interactions (LPI) and transport of X-rays from the cavity wall to the capsule, which involves the solving of complex kinetic 
and hydrodynamic equations implemented in the codes such as LASNEX [4]. In practice, simple mathematical models such as 
view-factor codes are usually employed to compute the radiation flux distributed on the capsule, especially for the preliminary 
design and optimization of thermonuclear target structure and shape [4], [5]. As described in [2], the radiation flux is usually 
obtained by solving a non-linear view-factor based equation model, which can be solved by utilizing Newton-Rapshon [6], 
Jacobi iteration [7], Cholesky decomposition [8], or Preconditioned Conjugate Gradient method [9]. However, in order to 
improve the evaluation accuracy, the size of discrete element is usually very small, the number of elements or equations will 
increase significantly, which may lead to much time required to solve large scale nonlinear equations. The computation time 
may be unacceptable for researchers. Therefore, a new efficient computation approach is essentially required to significantly 
improve the efficiency of radiation symmetry evaluation and optimization. 

Compressed sensing (CS), proposed by Donoho and Candes [10], [11], is a new method to reconstruct signals from 
significantly fewer measurements than that required by traditional methods, which has attracted considerable attention and 
achieves successful applications such as Medical imaging (MI) [12], Analog-digital Conversion [13], Computational Biology 
[14], and Computer Graphics [15]. In the field of radiation symmetry evaluation of ICF, compressed analysis approach has been 
applied to efficiently obtain the radiation flux distribution on the capsule [16]. However, the radiation flux computation model is 
simplified, and the accuracy of the radiation asymmetry evaluation is limited. The non-linear Time Dependent Energy Balance 
Model (TDEBM) presented in [17] is often used to compute the radiation flux. Therefore, compressive sensing approach for 
non-linear equation solving approach, such as Iterative Hard Thresholding (IHT) [18], can be applied. However, the step size 
is often taken as a constant, which may lead to more times of iteration convergence. Normalized Iterative Hard Thresholding 
(NIHT) [19] is further proposed, in which, fixed step length is replaced with a descending factor to accelerate the convergence 
rate. Nevertheless, the gradient based search direction in NIHT, may lead to the convergence is very slow when it approaches 
the minimum. Conjugate Gradient Iterative Hard Thresholding (CGIHT) is proposed in [20], to replace the search direction with 
the conjugate gradient direction, which can benefit accelerate the convergence. However, the search direction of CGIHT may not 
be conjugate, which may lead to more times of iterations being required. As one of the iterative greedy algorithms, Subspace 
Pursuit (SP) algorithm has attracted much attention for its backtracking idea in solving compressed sensing reconstruction 
problems, which requires only a small number of iterations [21]. However, SP algorithm can not be directly used to solve 
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Fig. 1: Scheme of the cylinder-to-sphere model and the energy exchange process for an element 


nonlinear problems. In this paper, we absorb the advantages of CGIHT algorithm and SP algorithm, and propose a Conjugate 
Gradient Subspace Thresholding Pursuit (CGSTP) algorithm to largely reduce the TDEBM and efficiently obtain the radiation 
flux. The core ideas include: leftmargin=8mm,labelsep=4mm 

(1) The spherical harmonics, annular Zernike polynomials and Legendre-Fourier polynomials are employed to sparsely 
represent the radiation flux on the capsule and spherical cavity, which enable only a few equations are required to 
recover the radiation flux with high accuracy, and significantly reduce the radiation computation model. 

(2) Over the three sets of polynomials domain, such reduced nonlinear equation model with sparse coefficients is formulated 
to enable compressed sensing algorithms, such as CGIHT, be used to efficiently obtain the radiation flux. 

(3) Conjugate Gradient Subspace Thresholding Pursuit algorithm is then presented to enable the adjacent search direction 
be conjugate, which facilitate the convergence of iteration, rapidly obtain the solution, and then efficiently evaluate the 
radiation symmetry for ICF experiments design. 

The outline of this paper as follows. Section 2 introduces the nonlinear view-factor based equation model in a cylindrical 
cavity model. Radiation sparse representation model constructed with three orthogonal bases is described in section 3. Then 
section 4 gives a description of the proposed reconstruction algorithm, i.e. Conjugate Gradient Subspace Thresholding Pursuit 
(CGSTP) algorithm. At last, section 6 concludes following experiments validation in section 5. 


II. SPARSE REPRESENTATION OF RADIATION FLUX 
A. Radiation flux computation model based on view-factor 


As shown in Figure 1 , eight laser beams are injected into the cylindrical cavity through two entrance holes. By interacting 
with the high-Z material on the inner wall of the cavity, the laser incident on the cavity is converted into X-ray, the X-ray 
irradiates evenly and drives the centrally located capsule implosion. With generated discrete mesh elements, the radiation flux 
on the wall of the cavity and surface of the capsule can be computed by using the view-factor function described in [16], i.e. 
the TDEBM. For a discrete mesh element, according to the energy conservation law, the energy received from other elements 
equals the energy emits to other elements. It can be formulated as: 


N 
A(E; + X. Fj:&GB;) = Bi, (1) 
j=l 
where JN is the total number of mesh elements, F; is the radiation flux comes directly from primary source converted from 
laser on the i-th mesh element which have a centroid p; and normal n;, can be computed as F; = N SOF, yCj, here, S9 
is the radiation flux converted directly from laser beams of the j-th element. Ç; is the area of the j-th mesh element with a 
centroid p; and normal n,. Fj; is the view-factor between the j-th mesh element and the i-th mesh element, which can be 
computed by 


4 
Py = E [i= (pj = pò]: Drs: (Pi = vy) / Ipi all D 


and À; is the albedo of the i-th mesh element, which is relating to the wall material, radiation time t and the emit energy B;, 
which can be given as 


1 
Mi = =] , 
14u87 BP tae 
where v , a and 8 are the material parameters in which v = 4.87, a = 8/13, 8 = 16/13. B; and B; are the radiation 
flux of the i-th mesh element and the j-th mesh element respectively. Let E = [F, E2,- +, ET , S = [99, 99,- -, so)", 


(3) 
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F = (Fy), ¢ = [Gai @,** Cul B = [Bi, Ba,:- By]; C =v teb, where T indicates the matrix transpose. Let 
V = (Fji¢;) be an intermediate variable, then the matrix form of E; = ye SP Fy; can be written as E = VS and 
Equation (3) can be written in a compact form as 


(I-V)B+CB°) = E (4) 


where I denotes an unit matrix. In Equation (4), V and E are known, B is the unique unknown variable. ° denote the 
hadamard power which means for any vector V and constant a, V°° = [V}, VS, =+, vay holds. Since 6 do not equal one, 
Equation (4) is nonlinear. Such nonlinear equations can be solved by utilizing some conventional iterative methods such as 
Newton-Raphson method. However, these methods is time consuming even unusable when the number of the discrete elements 
reaches 107 or more. 


B. Compressed Sensing Method 


Compressed Sensing (CS) also called Compressed Sampling, is an emerging technique, which has shown that sparse or 
compressible signal can be reconstructed from far fewer measurements than what is considered by Nyquist theorem [10]. 
Compressed sensing theory usually consists of sparse representation, measurement matrix, and sparse coefficient reconstruction. 

1) Sparse representation 


Consider a N x 1 signal f = [fi, fo,--- yin] which can be expressed as a linear combination of a complete basis W: 
N 
f=) abi = We (5) 
i=0 
where c = [c1, C2, ++, en], W = [y1, Ho,---, Wy]. If there are K non-zero elements in c, then c is called K-sparse and we 


say f is sparse in the © domain. 

2) Measurement matrix 

M projections of the signal f constitute a M x N(M < N) matrix ® called measurement matrix, then the sampling 
process can be described as: 


y= f =@®Vc= Ac (6) 


where A = ®W is known as measurement matrix. It is clear that the process of recovering c from y is an ill-posed problem 
because M < N and seems impossible to solve such underdetermined linear equation. However, if the sparse coefficients c 
is K-sparse, and ® satisfy the Restricted Isometry Property (RIP) criterion [22], then c can be accurately recovering with a 
high probability from M measurements. The RIP criterion can be expressed as 


(1 —dx)llellz < Iel] < (1 + ôx)llell. (7) 


The restricted isometry constant g (0 < g < 1) is defined as the smallest constant for which this property holds for all 
K-sparse vectors c. 

3) Sparse coefficients reconstruction algorithm 

As can be seen from above, the sparse coefficient vector c can be recovered with high probability from y by solving the 
following ọ-norm optimization: 


arg min||y — Acl|, s.t. ||e||>=K (8) 


where £o-norm, ||x||o, denote the number of non-zero entries in the argument and ||x||2 indicates its Euclidean norm. Some 
iterative greedy algorithms have been developed to solve such combinatorial optimization problem, such as Iterative hard 
thresholding (IHT) [18], Subspace Pursuit (SP) [21], and a series of algorithms derived from them. 


III. SPARSE REPRESENTATION OF THE RADIATION FLUX 


In the cylindrical cavity model of ICF, as shown in Figure 2, the radiation region located in the cavity consists of three 
parts, i.e., spherical surface of the capsule, the bottom and the top surface and the side surface of the cylindrical cavity. Here, 
B., B, and B, are used to represent the radiation flux of theses three parts respectively and the total radiation flux can be 
expressed as B = Bs + Bu + By. Since the X-ray converted from incident laser are uniformly and continuously distributed on 
the radiation regions, the radiation flux could be represented by linear combinations of some complete orthogonal polynomials. 
In this section, three orthogonal bases which are constructed by three sets of orthogonal polynomials are used to sparsely 
represent the radiation flux. 
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Fig. 2: Division of radiation surfaces in cylinder model in cartesian coordinate 


A. Sparse representation of the radiation flux on capsule 
Let (a ) py? be a spherical coordinate of the 7-th element, then the radiation Bs (08, oS) can be represented by a 
linear combinations of the spherical harmonics as follows [16]: 


B, (9, 0?) = aes CO (9) 


ns=0 


where the subscript ns is a polynomial ordering index starting with ns = 0, cn, is its coefficient and Y (0s, ps) is the spherical 
harmonics which can be represented as 


ba kepka cosk,yp, if kg >0 = 
Vine (9s, Ps) = Nav, Pg, (cos Os) { ike af heb o (—mg < kg < ms) (10) 
where NES = ,/(@ms+U(ms—ks)) is the normalization constant and PFs (cos 8s) is the Associated Legendre polynomial 
ms (4x(ms+ks)!) Ms s 8 poty : 


As shown in Figure 2, let AO, and Ay, be the minimal angles of subdividing the capsule over the domain 2(6,,y,) = 
(0, z] x [0, 27). Equation (10) can be transformed into an approximate form as 


Ns +00 
(9s, Ys) D 5 Cns Yn, (09, pP) ws 3 D Cn, Yn, (0% s ), of) =Ycf, (11) 
j=l ns=0 j=1 ns=0 
where Lg is the number of spherical harmonics expansion terms, c* = [c1, C2,- Chg]. is coefficient vector and Y is 


spherical harmonics basis which can be expressed as 


rie) deh) a 
Yo(6s z Ps ) Y (eF o ) TRS YL (of ys ) 
YNsxLs = . . Dr (12) 
YONO pres) Ce (Nee) siv Yr (0ND e) 
where Nọ, = Ao and No, = Ao denote the number of elements divided along the longitude and latitude respectively, and 


Ns is the total number of discrete elements on the capsule which can be calculated by Ns = Ng x Nọ. The 3D pictures of 
the first 16 terms of the spherical harmonics are shown in Figure 3. 


B. Annular Zernike Polynomials 


The bottom and the top surface of the cylindrical cavity is annular because it has two laser entrance holes. As described in 
[23], Zernike annular polynomials, defined on an annular region, is obtained from the Zernike circle polynomials [23] by using 
the Gram-Schmidt orthogonalization process. The Zernike annular polynomials are widely used to mathematically describe 
wavefront aberrations of optical systems due to it is orthogonal over an annular domain, and it can uniquely describe any 
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Fig. 3: 3D rendering of the first 16 terms of spherical harmonic polynomials 


continuous function with the same definition domain [22]. Let (ri! ) l Yu >) be a polar coordinate of the i-th element of the 


bottom surface, the annular Zernike polynomials have the following forms in a polar coordinate: 


ae (r, f) cos kupu for ky > 0 


E E ee 
Ur (r, 7; Pu) -{ RE, F) sin kupu for k, <0 a 


where r(0 < r < 1) and (f < r < 1) are the outer and inner semi-diameter of the annular region respectively, 
Pu (0 < yu < 27) is the circumferential angle, lu > 0, ku > 0 are integers, lu — |ku| > 0 and even, R“ (r, T) is the 
radial components of the polynomials. The orthonormality of the annular Zernike polynomials can be expressed as 


1 2m 1 2m 
os Fs py UR (r, searten | |  rrdon = Fuh (14) 


where 01,1", (Okun, ) stands for the Kronecker delta, equal to one if lz = I’; (kz = k’,) and to zero otherwise. when 7 = 0, 
the annular polynomials reduce to the Zernike circle polynomials. Let l, = 27 + ky, the radial polynomials expression of the 
annular Zernike polynomials can be obtained by the recurrence formula: [23] when ką = 0 


2 _ «2\ 1/2 2 _ p2 
RY = R(T ) -p PE A) (15) 


1—72 1—72 


where P;(-) is the Legendre polynomials. When k, > 0, the radial polynomials are given by the recurrence relationship 


kaa 1/2 
=f 
Ree nr= rii QE" (2) , 16 
Dik, (F) ao Q; (r°) (16) 


where Oy (T) is a set of orthogonal polynomials obtained by orthogonalizing the sequence 1, 7, +- ,7/ over the interval (F, 1) 
with a weight function 7%», H ~ are the normalization constants. the recurrence formula of the QF" (T) and H fu as follows: 


_ 2(2j+2ku—1) H E RETTORE (r?) 


Qk 2 — i = = ~ (17) 
J ( ) (j + ku) (1 — f2) OF 10) 7 He T 
; Kyi 
Hře — 2(2j + 2ku = 1) O51 (0) Hře! (18) 
I O (Jkl i) QETO) 
j 
especially, when ku = 0, 
Q90?) = RB (nF): H = RE (19 
j^) = Ry, T); H; = a7 +1)’ 


When k,, < 0, the radial polynomials can be obtained by the relationship: 
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Fig. 4: 3D rendering of the first 15 terms of Annular Zernike polynomials 


Bs ea i (ee). (20) 
When ku = lu 
ku 
RE (r,ř) = : z z 21) 
oe 


In practice, the outer diameter of the bottom surface is not equal to 1, so we require to map the reality interval [router, rinner] to 
the interval [7, 1]. Similar to the radiation flux representation on the capsule, the radiation flux Bu (re pe) can be expressed 


by a linear combinations of the annular Zernike polynomials as 


Ny +00 


ne)=S Y Su Un, (r°, pP) a SS ay mile, #0) = Uc", (22) 


i=l ny=0 g=] ny=0 


where nu is a polynomial ordering index starting with ną = 0, U is the annular Zernike polynomials with matrix form which 
presented as 


U(r, pl?) Uir O # pu) U(r, # pW) 
U(r, pu) ae oe) 4 Ur (r Fpl) 
Un, xt, = . . . . (23) 
Uo(r (N,) „F; ok Neu)) Uy (rN), ži oe) a Ur, (rN F; yN) 


where the subscript N, and N,, are the number of elements divided along the radial direction and rotation direction 
respectively. N, is the total number of elements of one of the end face of the cylinder cavity which can be calculated by 
Nu = N, X Ng,,. Ly is the number of the Zernike annular polynomials expansion terms which can be obtained by the formula 
Ly, = (Nu + 2) (nu + 1) /2. The 3D pictures of the first 15 terms of the annular Zernike polynomials are shown in Figure 4. 


C. Legendre Fourier polynomials 


Since the side of the cylindrical cavity is a cylindrical surface, it is no doubt that cylindrical coordinates (p, pw, z) seem 
to be the optimal choice for describing Bw, where p, Pw and Z denote radius, azimuth angle and height of the cylindrical 
coordinate. Obviously, Bw is only related to azimuth angle and height and not to radius. 

According to the construction principle of multivariable orthogonal polynomials, two sets of one-dimensional polynomials 
which are defined in the azimuth direction and height direction respectively, can be employed to construct a two-dimensional 
orthogonal polynomials for describing Bw. Since the cylinder surface profile along the azimuthal direction is closed and 
periodic, a Fourier series is suitable for characterizing the azimuthal coordinate, which can be expressed as 
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TABLE I: The first 9 expansion terms of the LF polynomials sequence 


Terms Order Polynomials Expression 
Wo 0 LoFo 1 
Wi 1 LoF_1 sinkwpw 
W2 1 LoFi cos kwpw 
W3 1 Lı Fo 5 
Wa 2 LoF_2 sin 2kw pw 
Ws 2 Lo F2 cos 2kw Pw 
We 2 Lı F—ı z sin kwpw 
W7 2 Li Fy z cos kwpw 
Ws 2 Lo Fo (32? — 1)/2 
cos k kw > 0 
Fru (Pu) = So ew Sw (24) 
snkypy ky <0 


where ku = 0,+1,+2,---. 

In height direction, coordinate z is restricted in the interval [—h,h], where h is half the height of the cylindrical cavity. 
A nature choice for z coordinate is a Fourier series [25]. Similar to the Associated Legendre polynomials used for spherical 
harmonics, some one-dimensional-polynomials such as Legendre polynomials and Chebyshev polynomials are suitable for 
the z coordinate. Since Legendre polynomials are orthogonal across the interval [—1, 1], which are conveniently mapped to 
the interval [—h, h], here we consider use Legendre polynomials for characterizing the z coordinate. Then a two-dimensional 
polynomials can be defined as the tensor products of Fourier series and Legendre polynomials in the variables z and Yw: 


Wry (2, Pw) = Pry (z) Few (Pw), (25) 


where P,,,(z) is Legendre polynomials with lẹ degree, Wn, (z, Pw) termed Legendre Fourier (LF) polynomials defined in the 

domain Q = {(z, pw) € [0,27] x [-1,1]}, nw denotes the ordering index of the LF polynomials starting with ną = 0. An 

integer nw represents the degree of the LF polynomials which can be calculated as ny, = lw +|kw|, where the operator |- | 

means taking the absolute value. It is clearly that the LF polynomials are separable in the cylindrical coordinates z and yw. 
Bw (z, Pw) can be approximately presented by the LF polynomials as 


Nw +00 Nw Lw 
Bu(z, Pw) = > > Cnu Wna (2®, pO) x 5 5 Cnu Wna (2®, 02) = Wc”, (26) 
i=l nyw=0 j=1 ny =0 


where cn, is the expansion coefficient of the LF polynomials.W denotes the LF orthogonal basis which can be written as 


Wo(2, 09) Me, pP) o Wee, pP) 
Wo(2, pP) Me, pP) Wr, (2, pA) 
WN, xL = . . , f j (27) 
Wole, pew) Wr (2%), pee) sea Wr, (202, pNew?) 


where the subscript N, and N,, are the number of elements divided along the radial direction and rotation direction 
respectively. N,, is the total number of elements of the side surface of the cylindrical cavity which can be calculated by 
Nw = Nz x No,- Lw is the number of the LT polynomials expansion terms. The top 9 expansion terms of the LT polynomials 
are given in Table I. The 3D pictures of the first 16 terms of the LF polynomials are shown in Figure 5. 


D. Sparse Representation Model of Radiation Flux 


From Subsection 3.1-3.3, the total radiation flux in the cavity can be expressed in the form of sum: 


B = B; + By + By = Yn,x1,c° + Un, x1,c" + WN, xL C”, (28) 
or matrix form: 
Yn, xL j 
B= UN,„x Lau x c” = We. (29) 
WN,,xL g" 


Equation (29) can be written as a functional form 
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Fig. 5: 3D rendering of the first 16 terms of Legendre-Fourier polynomials 


f(c) = (I1-V)(We) + C(We)ol/? — E (30) 


which has a unique unknown variable, i.e., sparse coefficient c. 
In this section, three sets of orthogonal polynomials are used to represent the radiation flux in a cylindrical cavity, and the 
radiation energy balance equation is transformed into an equation about the sparse coefficients. 


IV. RADIATION FLUX RECONSTRUCTION UNDER COMPRESSIVE SENSING FRAMEWORK 
A. Compressed Observation Model of Radiation Flux 


Compressed Sensing is usually works on linear system. First-order Taylor expansion around a point approximation is usually 
applied to linearization in nonlinear CS that is similar to the case in Equation (30). Equation (30) can be expanded at point 
c* as 


fle) fle") + fer(e— c), (31) 


where fe» is a linear operator such as Jacobian matrix which can be expressed as 


fe = (I- V)® + diag 30e- , (32) 


where diag(e) denotes a diagonal matrix whose elements on the diagonal are composed of e. Then Equation (32) can be 
rewritten as: 


fec = fec* — f(c*). (33) 


As the experimental results shown in Table 2, in the case of smaller grid, the computation is very large for solving view- 
factor matrix via Newton-Raphson method. In accordance to the CS theory, a measurement matrix of M x N(M « N) can 
be used for obtaining linear and non-adaptive measurements from Equation (33), which can be expressed as 


® foc = P(foc* — f(c*)). (34) 


The commonly used measurement matrices have stochastic Gaussian matrices, random Bernoulli matrices, partial orthogonal 
matrices, structured random matrices and cyclic matrices and sparse random matrices [11]. Since the basis W is orthogonal, 
from which the randomly sampled elements must meet the orthogonal, and also to meet the Restricted Isometry Property (RIP), 
the sampling points above the elements can be directly selected. This means only needs to calculation the sampling elements 
composed of measurement matrix. In practical, the measurement matrix is a similar diagonal matrix consisting of the sampling 
of row index of fe. 

Let A= ® fe, y = P (f.«c* — f (c*)), Equation (34) can be written as 


Ac= y. (35) 
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Since M < N, in general, recovery of c is a combination problem which is known to be NP-hard. Fortunately, under stricter 
conditions on A (i.e., Restricted Isometry Property), a range of different algorithms can be used to recover c from roughly 
O(K log(N)) measurements [10]. Based on the results of Section 3, we can roughly estimate the number of samples for 
capsule, the bottom and the top surface of the cavity and the side of the cavity, respectively. For example, for a model with a 
total elements number of 9776 (where the the number of elements for capsule, the bottom and top surface of the cavity, and the 
side of the cavity are 2592, 2016, 2016 and 3152, respectively), the number of samples for capsule is 30 x log(2592) = 102.4, 
the bottom and top surface have the same samples of 35 x log(2016) = 115.6, and the number of samples on the side of 
cavity is 100 x log(3152) = 349.9, and the total number of samples can be estimated as 150 + 150 x 2 + 400 = 850. The 
number of samples in other models is shown in Table III. 


B. Algorithms for solving model 


Sparse coefficient c can be recovered by optimizing the following sparse constraint problems 


é=argmin||y—Ac\|, s.t. ||el|,< K (36) 


some greedy iterative algorithms such as hard thresholding algorithms and orthogonal matching algorithms can be used for 
solving Equation (36). In this section, we proposed a new iterative greedy algorithm for solving the radiation flux compressed 
observation model based on two typical CS reconstruction algorithms. 

1) Conjugate Gradient Iterative Hard Thresholding Algorithm 

Iterative hard thresholding (HT) algorithm is usually used for CS reconstruction due to its simple operation and low 
computational complexity, and it’s iterative formula as 


Cn41 = Hx[en + pA'(y — Aen), (37) 


where Hx (e) is a nonlinear operator that sets all but the absolute value of largest K elements of c to zero, u is the step size 
and usually taken as a constant in IHT algorithm. The convergence of this algorithm was proven in [dd] under the condition 
that || A||2 < 1, in which case, the IHT algorithm converges to a local minimum. Since step size and the scale of sensing matrix 
have a great influence on the convergence of IHT [18], Normalized Iterative Hard Thresholding algorithm (NIHT) was proposed 
in to improve IHT algorithm [19]. The restriction of the snesing matrix A in NIHT is replaced by a descending iteration factor, 
which greatly improves the convergence and stability. However, the search direction in NIHT algorithm is gradient direction, 
it may result slow convergence near the minimum. Conjugate Gradient Iterative Hard Thresholding (CGIHT) algorithm was 
proposed in [20], in which, conjugate gradient direction was adopted as the search direction to speed up the convergence. 
The CGIHT algorithm can be stated as Algorithm 1, in which, Ø denotes an empty set and d denotes search direction, ~y is 
an orthogonalization weight, the support of c denoted by T,, and Ar, indicates a submatrix of A. By allowing the search 
direction change after each low complexity iteration, CGIHT is able to quickly identify the correct descending direction while 
using the computational advantages of conjugate gradient method when the support set remains stable. 


Algorithm 1 [20] Conjugate Gradient Iterative Hard Thresholding restarted for Compressive Sensing 
Input: A, y, K 
Initialization: cy = 0, T_, = 0, d_; = 0, To = supp(Hx(A7y)). 


for each iteration n > 0 do 
1) gn = AT(y— Acn) (compute the gradient direction) 
2) if TT, #fTn-ı then 
Xn =0 
else 
Xn = ge Ae (compute orthogonalization weight) 
end 
3) dn = 9Gn+Xndn—-1 (compute conjugate gradient direction) 
4) Qn = (Gnv,,9nt,)/(Ar, dnT,, Ar, dnr,) (compute step size) 
5) Cn+1 = Hg(£n F andn), Tn+1 = supp(€n41) 


until the stopping criteria is met 


end for 


Although the CGIHT algorithm can accelerate convergence, the Jacobian matrix needs updating at each iteration in non-linear 
problems, which may result in search directions that do not satisfy the conjugate condition, resulting in slower convergence. 
In the process of solving nonlinear radiation flux calculation model, it is time-consuming to update the Jacobian matrix, and 
we should to devise means to reduce the number of iterations. 
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2) Subspace Pursuit Algorithm 

Orthogonal matching pursuit (OMP) algorithm is wildly used for CS reconstruction owing to its simple implementation 
and lower computational complexity [21]. In each iteration, OMP selects the column of the measurement matrix which is 
most strongly correlated with the current residuals and adds this column into the set of selected columns, then the residual is 
updated by projecting the measurement vector y onto the linear subspace spanned by the selected columns. The reconstruction 
complexity of the standard OMP is roughly O(K MN) because it always needs s iterations. Unlike the OMP algorithm, at 
each iteration, subspace pursuit algorithm selected s candidates based on the correlation values between the columns of A and 
the measurement vector y, and then a back-tracking operation of the candidate set in the previous iteration is implemented. 
Compared to the OMP algorithm, the computational complexity of the SP algorithm can be further reduced to O(M N log(K)) 
when the nonzero entries of the sparse signal decays slowly. The SP algorithm can be described as follows. 


Algorithm 2 [21] Subspace Pursuit algorithm 
Input: A, y, K, 
Initialization: To = {K indices of the largest magnitude entries in the vector ATy}, 
Co = argmin,||y — Anella ro = y — Anco. 


When n > 1 Do 


1) T = Ta- U {indices of K largest magnitude entries of ATr,_1} 
2) ĉn = {arg min ||y — Aal|,,supp(x) € T} 

3) Tha = supp(Hx (én) 

4) cn = {argmin||y — Ax||,,supp(x) € Tr} 

5) Tn =y— Ar, en 


until the stopping criteria is met 


Output: cn, Tn 


3) Conjugate Gradient Subspace Thresholding Pursuit Algorithm 

Inspired by CGIHT algorithm and SP algorithm, a new greedy iterative algorithm named conjugate gradient subspace 
thresholding pursuit (CGSTP) is proposed for reconstructing radiation flux. This algorithm is based on SP algorithm, drawing 
on CGIHT algorithm and using conjugate gradient as the search direction to accelerate convergence. The pseudo-code of 
CGSTP algorithm is shown in Algorithm 3. 


Algorithm 3 Conjugate Gradient Subspace Thresholding Pursuit (CGSTP) 
Input: A, y, K 
Initialization: co = 0, T_1 = 0, d—ı = 0, To = supp(Hx(A7y)). 
for each iteration n > 0 do 
1) gn = AT(y— Acn) (compute the gradient direction) 
2) if T, #Tn-ı then 


Xn =0 
else 


= (AgnTy:Adn—1Tp) 
Xn = (Ady —17,,,Adn—1T,, ) 
end 


3) dy = Qnt+Xndn—-1 (compute conjugate gradient direction) 
4) Hn = (QnT,, , gnr, )/ (AT, dnt, , Ar, dnr, ) (compute step size) 


5) Thi = supp(Hx (en + Lndn)) UT). 
6) én41 = {argmin||y — Aal|,,supp(x) € Trai}, Troi = supp(én41) 


(compute orthogonalization weight) 


T) Cn41 = {arg min ||y — Ag||,,supp(x) € Tr+41} 
until the stopping criteria is met 


end for 


The CGHTP algorithm is a simple combination of the CGIHT algorithm and the SP algorithm, and absorbs the advantages 
of both algorithms. SP algorithm introduces backtracking operation in the iteration process. In each iteration, the support set of 
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the previous iteration is reconsidered, and the new support set is updated by the least square method. For compressible signals, 
SP algorithm can find the best s atoms in log(K) iteration. However, SP algorithm has strict requirements for RIP conditions. 

The CGIHT algorithm optimizes the search direction, and the conjugate direction is used as one of its optimization directions, 
which has a superlinear convergence rate associated with (vK — 1)/(./« + 1) [26], where « is the condition number of matrix 
A. However, it does not make further accurate estimation of the sparse coefficient after updating the support set, and the 
CGIHT algorithm only shows its advantage in convergence speed when the correlation of the columns in A is very small. 

The CGIHT algorithm uses the negative gradient direction and the conjugate gradient direction as the search direction. If 
the sensing matrix remains unchanged in all steps, only a small number of iterations are needed to converge. In fact, since the 
incomplete orthogonality between the columns of the sensing matrix, the support set selected in each iteration may be different, 
and the sub-matrix formed by the columns of A with the support set as the index is different in each iteration, which will 
lead to the search directions obtained in the previous iteration and the current iteration are not conjugate, thus the superlinear 
convergence speed of the algorithm can not be ensured. What’s more, in the nonlinear problem, because the CGIHT algorithm 
does not accurately estimate the sparse coefficients in each iteration, the solution process requires more iterations. 

The CGSTP algorithm absorbs the advantages of SP algorithm and CGIHT algorithm in support set update strategy and 
search direction optimization respectively. CGSTP provides two alternative search directions in the iterative process, including 
the gradient descent direction and the conjugate gradient direction, and adaptively selects the optimal direction according to 
the criterion of whether the support set is the same in the previous iteration and the current iteration. In theory, if the sensing 
matrix in all iterations remains unchanged, CGSTP can achieve superlinear convergence rate [26]. However, in solving the 
non-linear compressed observation model of radiation flux, the Jacobian matrix changes in each iteration, which results in 
the change of the sensing matrix, and the acceleration effect is not obvious. Fortunately, since CGSTP uses a backtracking 
operation similar to SP algorithm to accurately estimate the sparse coefficients after each updating of the support set, resulting 
in a small change in the Jacobian matrix obtained from the previous iteration and the current iteration, which increases the 
probability that the sensing matrices in different iterations are the same. 


C. Overview of radiation flux computation process with CGSTP 


According to the above analysis, solving the non-linear radiation energy balance equation in the Compressed Sensing 
framework can be summarized into three modules, namely, the sparse representation module of radiation flux, the compressive 
sampling module and the reconstruction module. Figure 6 summarizes the basic contents of each module. 

As can be seen from Figure 6, in the sparse representation module, all the radiation regions of the cylindrical cavity model are 
first divided into discrete elements according to a given precision, and parameter data such as coordinate values, normal vectors, 
and spatial angles are obtained. Then, the sparse basis composed of orthogonal polynomials of different expansion orders is 
calculated according to the elements parameters. The representation error of radiation flux over the sparse basis is analyzed 
to determine the number of expansion terms of polynomials, and the sparsity of radiation flux in the sparse transformation 
domain is determined by simulation experiments. Finally, a sparse representation model of radiation energy for a small number 
of sparse coefficients is constructed. 

In the compressive sampling module, firstly, all discrete mesh elements are indexed, and M elements are collected by random 
uniform sampling method (Latin Hypercube Sampling (LHS) method used in this paper). Then, the sparse representation model 
of radiation flux is transformed into a compressed observation model by using the equation corresponding to M sampling 
indexes, so as to realize the compression and dimensionality reduction of the model. At the same time, the occlusion condition 
of each element in the cylindrical cavity is analyzed, and the partial view-factor matrix is calculated according to the sampling 
index. 

In the reconstruction module, the linear approximation system is first constructed by using the first-order Taylor expansion at 
the current iteration point to construct a linear approximation system. Then, the partial Jacobian matrix is calculated according 
to the sampling index, and the sensing matrix is calculated according to the Jacobian matrix, The minimization problem based 
on o norm is established when the sparsity is known as s. Finally, the above Zọ norm minimization problem is solved by 
CGSTP algorithm. 


V. SIMULATION EXPERIMENTS AND DISCUSSION 
A. Experimental model and simulation environment 


In this section, two cylindrical cavity models for Shenguang II and Shenguang-III laser devices are used for simulation 
experiments to verify the performance of the proposed method. The cylindrical cavity model in Shenguang II is shown in 
Figure 7. The models of Shenguang-II and Shenguang-III are represented by ”S2” and ”S3” respectively. All experiments are 
tested in MATLAB R2015b [27], running on a Windows 10 machine with Intel 15-7500 CPU 3.4GHz and 8GB RAM. 

The elements parameters of the S2 model and the S3 model are shown in Table II, where S2-1, S2-2, S3-1, and S3-2 
represent models having different elements sizes and their total elements numbers are respectively 9776, 38952, 20736 and 
82944. 
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Fig. 7: A cylinder-to-sphere model on the S2 laser facility. (a) S2 laser facility, (b) CAD model of S2, (c) radiation energy 


computed on the capsule, (d) unfolding view of radiation energy on the capsule 
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The random uniform sampling method used in this section is the Latin Hypercube Sampling (LHS) method. The sampling 
number is estimated according to the formula slog(N). Taking the S21 model as an example, the sampling numbers of three 
radiation regions are calculated as 30 x log(2592) ~ 102.4 (take 150), 35 x log(2016) = 115.6 (take 150) and 100 x log(3152) = 
349.8 (take 400). The sampling rates of different models are shown in Table II. 

As shown in the Table II and Table III, we can see that (1) the number of discrete mesh elements significantly increases 
with the decrease of the size of mesh elements, and (2) the sampling rate for Compressed Sensing based approaches gradually 
decreases as the number of mesh elements increases. 


TABLE II: The size and number of discrete elements of four simulation models. ”H”, ”C”, ” Ar” and ”Ah” represent the 
minimum dimensions of the longitude direction of capsule, the latitudinal direction of capsule, the radial direction of cavity, 
and the axial direction of cavity, respectively. 


Model Capsule The end face of cavity The cylinder surface Total number 
2592 2016x2 3152 
Sal (H:5°, C:5°) (H:2.5°, Ar:15um) (H:5°, Ah:30um) 9776 
10368 8064x 2 12456 
52-2 (H:2.5°, C:2.5°) (H:1.25°, Ar:7.5um) (H:2.5°, Ah:15um) 3829% 
2592 4320 9504 
Sael (H:5°, C:5°) x 2(H:2.5°, Ar:15um) (H:5°, Ah:30um) 20736 
10368 17280x2 38016 
S32 (H:2.5°, C:2.5°) (H:1.25°, Ar:7.5um) (H:2.5°, Ah:15um) oer 


TABLE III: Number of samples of four simulation models 


Model Capsule The end face of The cylinder Total sampling Sampling rate 
cavity surface number 
S2-1 150 300 400 850 8.7% 
S2-2 150 300 450 900 2.3% 
S3-1 150 300 350 800 3.9% 
S3-2 150 400 400 950 1.1% 


B. Sparse representation error analysis of radiation flux 


Firstly, we verify the sparsity level of radiation flux Bs, Bu and B,, over three orthogonal bases Y, U and W, respectively. 
The expansion coefficients of three sets of polynomials and the representation errors (root mean squared error), between the 
original flux and the recovered radiation flux with ny, nu and nwu expansion order are shown in Figure 8(a) and (b) respectively. 
The relationship between the sparsity level of the coefficients and the representation error of the radiant energy flow is shown 
in Figure 9(a) and the magnitude of the sparse coefficients are presented in 9(b). 

As shown in Figure 8(a) and (b), with the increase of the number of expansion terms of orthogonal polynomials, the 
representation error decreases gradually. When the expansion terms of Spherical Harmonic polynomials, Annular Zernike 
polynomials and Legendre-Fourier polynomials reaches 400, 325 and 1225 respectively, the representation errors approach to 
zero and the curve remains stable. This means that only a small number of expansion terms are needed to accurately represent 
the complete radiation flux. From Figure 9 (a) we can see that the number of the absolute value of coefficients larger than 
107° over three sets of orthogonal bases are s1 ~ 30, s2 ~ 35 and s3 ~ 100 respectively. 

Since the top 35 order energy accounts for more than 99.9% of the total energy and the energy of more than 35th order 
is little. That is the analysis of radiation fluxs harmonic expansion situation under various irradiation conditions. The drive 
asymmetry must be less than 2% and the calculation accuracy should be setting at least an order for magnitude, equal to 0.1%. 


C. Solution results of radiation energy balance equations 


The resulting radiation flux on the capsule is computed with TDEBM and then compared with Newton-Raphson (NR), 
Preconditioned Conjugate Gradient (PCG), IHT, CGIHT and CGSTP on computation time, accuracy, iteration number for 
different kinds of mesh elements sizes. In all tests, the IHT, cgiht, and cgstp algorithms were run independently 20 times, and 
the average of 20 results was used as the test result. Then the computation time list for the four calculation models using 
different algorithms is shown in Table IV and V, and the iteration steps and accuracy of different algorithms are shown in 
Table VI and Table VII respectively. 

As shown in Table IV and Table V, we can see that (1) as the number of non-linear equations exponentially increases, 
the computation times thus increase sharply for traditional approaches, (2) less computation time is required for compressed 
sensing based approaches since the computation model is significantly reduced, only almost 8.7%, or 1.1% of equations are 
utilized to obtain the solution, and (3) since the search direction is optimized and the updating strategy of support set is 
improved, less time required for CGSTP than that of IHT and CGIHT to obtain the solution. Figure 10(a) shows the time 
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Fig. 8: The relation of representation errors and number of expansions over orthogonal polynomials. ”Z 2016”, ”L 3152” and 
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ratio between the CGSTP algorithm and the PCG algorithm for calculating the four models. Figure 11(b) shows the iteration 
time difference between the CGSTP algorithm and the IHT and CGIHT algorithms. It can be seen from Figure 10(a) that as 
the model scale increases, the time acceleration ratio of the CGSTP algorithm relative to the PCG method becomes larger and 
larger, and the maximum acceleration ratio can reach 80. This is because as the scale of the equation increases, the sampling 
rate becomes smaller and smaller, and the required calculation amount is less and less than the total calculation amount. Figure 
10 (b) shows that the difference between CGSTP and IHT and CGIHT is larger and larger with the increase of model size. 
This is because the radiation flux tends to be continuous with the increase of the number of discrete mesh elements, which 
reduces the number of updates of Jacobian matrix. This means that CGSTP has more advantages when the scale is larger. 
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TABLE IV: Running time of conventional method (unit: second),’—’ denotes out of memory 
Newton-Raphson method PCG method 
View-factor calculation Iteration Total Iteration Total 
Model ` ; z ; : 
time time time time time 
S2-1 4341.5 27.3 4368.8 21.4 4362.9 
S2-2 69045.8 2203.6 71249.4 313.2 69359.1 
S3-1 16253.1 1288.4 17541.5 96.8 16349.9 
S3-2 255843.2 — — 1572.6 257415.8 
TABLE V: Running time of Compressed Sensing method (unit: second). 
IHT CGIHT CGSTP 
View-tactor Sparse ba als Iteration Total Iteration Total Iteration Total 
Model calculation calculation : ; : ; : É 
: z time time time time time time 
time time 
S2-1 407.1 11.8 332.4 751.3 167.9 586.8 13:7 432.6 
S2-2 1518.4 40.2 1142.6 2701.2 683.1 2241.7 38.6 1597.2 
S3-1 658.6 Dont 671.8 1353.1 322.4 1003.7 20.7 705.0 
S3-2 3014.2 103.5 2248.3 5366.0 1131.2 4248.9 84.5 3202.2 
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Fig. 10: Running time of different algorithms. (a)The ratio of solving time between CGSTP algorithm and PCG algorithm, 
(b)ATS1 and ATS2 represent the difference of iteration time between CGSTP and CGIHT, CGSTP and IHT algorithm, 


respectively. 


The convergence performance of the above algorithms is shown in Table VI.It can be seen that (1) a small number of 
iterations is required for the NR, PCG and CGSTP approaches, whereas hundreds of iterations is required for the IHT, CGIHT 
approaches, (2) Even though less time of iterations are required for the NR and PCG approaches, the computation time for 
each iteration is much longer than that of IHT, CGIHT and CGSTP algorithms, and (3) as compared with the IHT and CGIHT 
algorithms, less iterations is required for the CGSTP approach since optimal search direction and step length is applied. Figure 
11 shows the relationship between IHT, CGIHT and CGSTP algorithm iteration times and reconstruction error. From Figure 
11, we can see that the error curve of CGSTP algorithm decreases very fast, and it only needs about six iterations to approach 
zero and keep stable, which means that the CGSTP algorithms has a faster convergence speed than the other algorithms. 


TABLE VI: Number of samples of four simulation models 


Model Newton-Raphson PCG IHT CGIHT CGSTP 
S2-1 3 10 1273 841 38 
S2-2 2 8 1124 729 35 
83-1 3 9 837 657 31 
$3-2 2 8 811 594 21 


The accuracy of the different algorithms is listed in Table VII. It can be seen from Table VII that (1) the NR algorithm and 
the PCG algorithm have the highest accuracy, and the root mean square error (RMSE) is lower than 1075, (2) the accuracy 
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Fig. 11: Iteration number and reconstruction error of different algorithms 


of CGSTP algorithm is higher than IHT algorithm and CGIHT algorithm in solving different models, and (3) the accuracy of 
all algorithms is improved with the increase of the number of mesh elements. The RMSE of CGSTP algorithm in calculating 
different models is less than 107, which satisfies the radiation symmetry analysis. 


TABLE VII: Reconstruction RMSE of different methods 


Model NR PCG IHT CGIHT CGSTP 
S2-1 3.14 x 10719 5.52 x 1076 3.71 x 1073 1.53 x 1073 7.84 x 1074 
S2-2 2.33 x 10719 3.91 x 1076 1.27 x 1073 9.84 x 107-4 6.20 x 1074 
S3-1 1.65 x 10719 4.13 x 1078 9.81 x 107-4 8.92 x 107-4 5.75 x 1074 
S3-2 1.08 x 10719 1.64 x 10-8 9.55 x 1074 8.51 x 1074 5.24 x 1074 


VI. CONCLUSIONS 


The radiation flux distribution symmetry evaluation is very important in ICF experiments, which involve a time consumption 
process in solving the nonlinear energy equilibrium computation model. In order to accelerate such equation solving process, 
an efficient radiation computation approach is presented. Firstly, we investigated the distribution characteristics of radiation 
flux on a cylindrical cavity model and three sets of orthogonal polynomials are employed to accurately and sparsely represent 
the radiation flux. Then, only a few mesh elements are sampled to formulate the sparse equation model, which enable the non- 
linear equation model largely reduced. Finally, a greedy iterative algorithm named CGSTP is presented to efficiently solve the 
non-linear radiation energy balance equations. Some experimental targets is utilized to validate the efficiency of the presented 
approach. The result show that the radiation flux can be efficiently obtained in almost 1/80 time of traditional approach. In 
future, we will continue to explore the sparse representation for recent free-form hohlraums, to enable the present approach 
be used for more targets design, analysis and optimization. 
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